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1.  Introduction 

Acetylcholinesterase  (AChE)  helps  to  regulate  some 
of  the  most  fundamental  human  biological  operations 
such  as  motor  function,  sleep,  attention,  memory  and 
emotions.  A wide  variety  of  AChE  inhibitors  exist,  with 
health  implications  that  range  from  fatal  to  benign. 
Strong  AChE  inhibitors  such  as  some  organophosphorus 
(OP)  nerve  agents  (some  of  which  can  phosphorylate 
AChE  irreversibly)  are  among  the  most  powerful 
neurotoxins  known  and  pose  severe  human  health  risks  if 
misused.  Weaker  ligands,  however,  have  been  broadly 
studied  as  possible  therapeutics  for  neural  dysfunction 
disorders  such  as  myasthenia  gravis  (MG),  and  various 
dementia  such  as  Alzheimer’s  disease  (AD).  Weaker 
inhibitors  have  also  been  employed  for  nerve  agent 
prophylaxis,  although  safety  concerns  continue  to 
surround  the  standard  prophylaxic  formulations, 
especially  pyridostigmine  bromide  (PB)^'',  suggesting  that 
the  criteria  for  their  approval  for  military  use  might  not  be 
acceptable  for  broader  civilian  medical  application.  In 
fact,  demonstrably  safe  and  effective  nerve  agent 
prophylaxes  remain  elusive,  despite  a wealth  of  AChE 
inhibition  data  arising  from  AD  and  MG  pharmaceutical 
research,  as  well  as  from  pesticide  toxicology  studies. 
Prophylaxis  design  is  inherently  more  challenging  than 
simple  inhibitor  design,  as  one  must  find  a balance 
between  an  interaction  strong  enough  to  prevent  adhesion 
of  other  ligands  but  weak  enough  to  permit  return  to 
normal  enzyme  function  in  a reasonable  time  frame. 
These  dynamic  and  competitive  factors  make  prophylaxis 
difficult  to  optimize  via  most  standard  in  vitro  assays  but 
lend  themselves  very  well  to  molecular  simulations. 

This  project  uses  a variety  of  modeling  techniques  to 
AChE  inhibitor  kinetics  and  dynamics  to  elucidate  and 
evaluate  prospective  new  prophylaxes  for  nerve  agent 


toxicity  threats.  Three  levels  of  theory  are  employed  to 
address  the  following  aspects  of  this  extremely  complex 
problem: 

Virtual  screening  of  inhibitor  formulation: 
Multiple  AChE  inhibitor  scaffolds  are  already  known, 
thus  ensuring  that  sizeable  virtual  libraries  of  potential 
prophylaxes  can  be  constructed  using  conventional  drug 
design  techniques.  The  computational  rigor  of 
prophylaxis  evaluation  requires  extensive  prescreening 
via  methods  such  as  Quantitative  Structure  Activity 
Relationships  (QSARs)  and  molecular  docking.  To 
accomplish  this  for  the  challenging  AChE  binding 
environment  (multiple  binding  subsites  and  ligand 
conformational  modes),  we  have  employed  three 
dimensional  (3D)  QSAR  type  correlations  to  formulate  a 
weighted  residue-based  scoring  methodology  specifically 
tailored  toward  accurately  evaluating  AChE  docking 
results  and  formulating  de  novo  inhibitors  based  on  the 
resulting  3D  receptor  interaction  field. 

Analysis  of  competitive  transport  kinetics: 
Fundamental  to  the  high  AChE  enzymatic  efficacy  is  a 
collaborative  network  of  receptor  subsites  that 
dynamically  guide  ACh  (or  inhibitors)  toward  the 
esteratic  subsite^^'.  Comparing  the  transport  properties  of 
different  ligands  and  their  competition  for  access  to  (and 
transit  past)  these  subsites  will  provide  insight  into  the 
requirements  for  effective  prophylaxis.  Such  analysis,  as 
well  as  quantification  of  transport  kinetics,  is  readily 
accomplished  through  MD  simulations  of  prophylaxes 
and  toxins  (alone  or  in  competition)  interacting  with  the 
solvated  receptor  environment. 

Prediction  of  adsorption  kinetics, 
thermodynamics,  and  physical  properties:  Predicting 
reaction  kinetics  and  thermodynamics  of  nerve  agents  and 
prospective  prophylaxes  requires  the  use  of  QC  methods. 
This  data  is  used  to  predict  bound  state  geometries  and 
transition  states  for  receptor-ligand  complexes,  and  to 
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employ  the  corresponding  complex  energies  to 
characterize  their  associated  thermodynamics  and 
chemical  kinetics.  Much  time  and  effort  has  been 
expended  in  delineating  the  precise  interactions  within  the 
active  site  and  neighboring  residues  upon  binding. 
Besides  the  obvious  intrinsic  worth  of  data  providing 
detailed  information  on  the  structural  features  which 
stabilize  binding,  unbinding,  and  related  side  reactions 
(such  as  aging),  these  results  provide  an  excellent  starting 
point  for  mutagenesis  studies  to  provide  further  insite  into 
agent  binding.  QC  studies  have  also  been  performed  to 
analyze  the  gas  phase  behavior  of  hydration  complexes  of 
nerve  agent  simulants,  as  a step  toward  understanding  the 
ambient  volatility  of  these  systems. 

2.  Computational  Details 

Virtual  Screening  of  Inhibitor  Formulation:  The 
AChE  crystal  structures  used  herein  were  obtained  from 
the  PDB  database  (http://www.rcsb.org).  Our  main 
docking  and  training  activities  have  been  focused  on  a 
human  AChE  (huAChE)  structure  (code#:lB41^^'),  but 
rely  on  ligand  binding  information  from  a torpedo 
California  AChE  (tcAChE)  structure  (code#:lEVE^‘'*)  that 
includes  a co-crystallized  E2020  inhibitor.  All  waters 
were  removed  from  the  structures.  The  huAChE  structure 
was  aligned  to  the  tcAChE  in  SYBYL^^’  to  ascertain  the 
orientation  of  the  ligand  E2020  relative  to  the  huAChE 
structure.  The  Root-Mean  Square  Derivation  (RMSD) 
between  the  two  is  only  0.85  A for  all  backbone  C^  atoms 
in  those  structures,  suggesting  good  alignment  and 
structural  similarity.  Hydrogen  atoms  were  added  (via 
SYBYL)  to  the  resulting  huAChE-E2020  complex.  The 
positions  of  these  new  protons  were  then  optimized  in 
MOE'®’  via  molecular  mechanics  using  the  MMFF94s 
force  field  (all  heavy  atoms  fixed)  to  avoid  bad  inter- 
atomic contacts.  The  position  of  E2020  was  then 
optimized  (all  receptor  atoms  fixed)  to  determine  a 
plausible  stable  conformational  structure  for  the  ligand  in 
the  receptor  environment.  In  both  of  the  above 
simulations,  MMFF94s  charges  were  used  to  account  for 
relevant  electrostatics.  The  Steepest  Descent 
minimization  algorithm  was  used  for  the  first  100  steps 
(unless  a RMS  gradient  of  less  than  100  kcal/(mol  A)  was 
first  achieved),  followed  by  200  steps  of  Conjugate 
Gradient  (unless  a RMS  gradient  of  less  than 
1 kcal/(mol  A)  was  attained),  and  finally  completed  by 
1000  steps  of  Truncated  Newton  (or  a RMS  gradient  of 
less  than  0.01  kcal/(mol  A)).  The  resulting  E2020 
structure  was  then  extracted  for  subsequent  docking 
calculations. 

A collection  of  69  compounds  with  IC50  data 
measured  with  human  AChE  assay^’*’^'  was  selected  for 
training  and  testing  the  scoring  function.  The  activity 


among  these  compounds  ranges  from  0.33  nM  to  30,000 
nM.  Docking  calculations  were  carried  out  with  the 
GOLD  program**®'.  A genetic  algorithm  was  used  in 
searching  the  binding  conformation  of  flexible  ligands, 
using  the  default  parameters  in  GOLD.  The  active  site  for 
the  huAChE  docking  calculations  was  constructed  from 
the  crystal  structure  by  retaining  all  residues  within  a 
radius  of  12  A relative  to  E2020  (but:  discarding  the 
original  ligand  itself).  A maximum  of' 20  poses  were 
computed  for  each  compound.  ^ose  docked 
conformations  were  saved  in  SDF  format  and  then 
imported  into  SYBYL  for  scoring  calculations  according 
to  the  FlexX  and  CSCORE  modules.  1 The  scoring 
methods  available  included  empirical  methods  such  as 
Chem  score,  FlexX  score,  and  G score,  ^d  knowledge- 
based  methods  such  as  PMF  score  and  Drug  score.  Multi- 
Linear  Regression  (MLR)  was  used  to  obtain  a consensus 
score  from  these  methods.  One  conformation  was 
selected  for  each  compound,  so  as  to  give  a good 
compromise  between  the  best  consensus  score  and  those 
with  the  closest  alignment  to  the  original  E2020  ligand. 
The  chosen  conformations  were  used  to  fit  an  interaction 
field  of  the  following  form:  | 

pIC5Q  = (1) 

' ^ I 

Where  c„  dj  are  fitted  coefficients,  Ef^  are 

electrostatic  and  VDW  interactions  arising  between  atoms 
in  the  residue  and  the  ligand.  In  this  expression,  all 
receptor  residues  within  10  A of  the  position  of  the 
original  E2020  ligand  (a  total  of  92  [residues)  were 

included  in  the  summation  over  '/,  and  Ef^ 

and  E''j^'^  calculated  via  an  SVL  script  written  for  the 

MOE  system.  The  statistic  analysis  was  performed  in 
Simca-P'"'  with  partial  least  square  regression  (PLS). 
Fifty-three  of  the  full  69  compounds  were  selected  as  our 
training  set,  and  the  other  16  compounds^  were  used  as  a 
test  set  for  validating  the  predictive  power  of  the  new 
scoring  function.  ' 

Competitive  Transport  Kinetics:  In  order  to  arrive 
at  starting  structures  for  simulations  | of  the  above 
collection  of  ligands  interacting  with  the  tcAChE 
receptor,  bond  distances  and  angles  for  each  ligand  were 
optimized  via  quantum  chemical  | HF/6-31G(d,p) 
calculations.  The  resulting  structure  was  'then  allowed  to 
torsionally  relax  within  the  receptor  environment  and  was 
generated  via  FlexX  molecular  docking  simulations.  In 
each  case,  the  ligand-receptor  complex  with  the  top  score 
(out  of  50  docking  runs  per  ligand)  as  determined  by  the 
Drug  Score  evaluation  scheme  (along  with  a visual  sanity 
check)  was  used  as  a starting  point  for  further  simulations. 
The  ligand — enzyme  complex  was  then  prptonated  via  the 
X-LEaP  program  distributed  within  the  AMBER  6 


software  suite^'^'.  For  all  calculations,  glutamate  and 
aspartate  monomers  were  modeled  in  their  anionic  form, 
lysines  were  described  as  cations,  N®  on  histidines  were 
protonated,  while  N®  were  left  in  unprotonated  sp^ 
hybridized  form.  All  other  residues  were  left  in  then- 
neutral  form,  except  (for  one  set  of  calculations)  the 
deprotonated  ser  200  anion.  Negative  charges  were 
balanced  with  Na+  ions  added  according  to  a coulomb 
potential  on  a 1.0  A grid.  In  each  case,  the  resulting 
complex  was  then  fully  solvated  within  an  aqueous  box 
defined  by  a solvent  buffer  of  7 A using  the  TIP3P  model. 
To  ensure  solvent  conservation,  three  dimensional 
periodicity  was  effected  via  an  Ewald  summation.  RESP 
atomic  charges  for  the  ligand  (and  the  anionic  serine 
residue)  were  derived  fi-om  Hartree-Fock/6-31C}(p,d) 
quantum  chemical  calculations  using  the  antechamber 
module  in  amber7.0.  Atomic  charges  for  all  other 
enzyme  atoms  corresponded  to  the  default  charges 
available  in  the  Amber  99  force  field.  All  subsequent 
classical  simulations  were  effected  with  the  AMBER  6 
program,  using  the  Amber  99  force  field.  The  solvated 
receptor-ligand  complex  was  minimized  via  standard 
molecular  mechanics  for  1000  steps  in  order  to  eliminate 
unrealistic  atomic  contacts.  The  system  was  then 
equilibrated  at  constant  pressure  fi-om  0-300  Kelvin  via  a 
20  ps  variable  temperature  MD  simulation  with  a time 
step  of  2 fs.  This  relaxed  the  ligand-receptor  complex  and 
its  solvent  shell,  permitting  the  box  to  swell  to  dimensions 
of  slightly  less  than  80  x 80  x 80  A^  Equilibration  also 
has  the  effect  of  randomizing  the  ligand  conformation  and 
position  within  the  receptor.  From  this  relaxed  and 
randomized  structure,  dynamic  analysis  on  the  system 
was  effected  via  a 1 ns  constant  volume  simulation  in 
which  the  system  was  warmed  fi-om  100-300  K.  In  each 
case,  the  entirety  of  the  trajectory,  except  approximately 
the  first  30  ps  of  the  simulation,  persisted  within  a 295- 
305  K temperature  window  that  permitted  reasonable 
isothermal  analysis. 

Prediction  of  Adsorption  Kinetics, 

Thermodynamics,  and  Physical  Properties:  All 

calculations  are  performed  with  the  Gaussian  suite  of 
codes^'^^.  Starting  structures  for  acetylcholinesterase  are 
obtained  fi-om  X-ray  crystallographic  data  and 

coordinates  for  the  residues  of  interest  are  extracted  using 
the  RASMOL  code.  The  largest  models  incorporate  the 
active  site  residues  (Ser200,  His440,  and  Glu327  in  the 
Torpedo  Califomica  numbering),  the  oxyanion  hole 
residues  (Glyll8,  Glyll9,  and  Ala201),  and  Glul99. 
Calculation  are  performed  at  the  using  density  functional 
theory  B3LYP  (unless  otherwise  noted)  with  the  3-2 IG 
basis  set.  These  results  are  then  used  as  the  initial 

structure  for  optimization,  fiequency  and  NMR 

calculations  performed  at  b31yp/6-3  lg(d,p),  unless 
otherwise  noted. 


3.  Results 

Virtual  Screening  of  Inhibitor  Formulation:  Our 
scoring  model  built  via  PLS  regression  over  interactions 
within  the  53  molecule  training  set  appears  to  be  of 
reasonable  quality,  with  a correlation  coefficient  of 
R^=0.89  and  a leave-one-out  cross-validation  correlation 
of  Q^=0.72.  In  using  the  scoring  function  to  evaluate  the 
activity  for  the  16  molecules  testing  set,  we  achieved 
good  predictivity,  a correlation  of  R^=0.69  between  the 
calculated  results  and  experimental  values. 

To  compare  the  precision  and  extensibility  of  our 
scoring  function,  we  contrasted  the  above  results  with 
predictions  made  using  several  commercially  available 
scoring  methods,  including  Chemscore,  Flex  score.  Drug 
score,  G score  and  PMF  score.  The  correlation  between 
the  experiment  and  any  single  score  method  is  poor.  The 
PMF  score  showed  the  best  correlation  but  was  still  rather 
poor  (R^=0.13  for  the  training  set).  All  of  the  other 
representations  gave  even  worse  correlations;  ChemScore 
= 0.07,  FlexX  = 0.05,  DrugScore  = 0.01,  G Score  = 0.004. 
Since  those  general  scoring  methods  were  not  trained  in 
this  AChE  inhibition  set,  it  is  not  completely  fair  to 
compare  them  directly  with  our  specially-tailored  energy 
decomposition  method.  Therefore,  we  built  a consensus 
score  by  training  (over  the  53  molecule  set)  a weighted 
sum  over  the  above  five  commercially  available  scoring 
functions  as  follows: 

pIC50=  0.05682  * Chemscore  -0.00499*Drugscore- 
0.03582*Flexscore 

-0.01232  ♦ Gscore  +0.01847  ♦PMFscore  +7.62820  (2) 

where  Chemscore,  Drugscore,  Flexscore,  Gscore  and 
PMFscore  refer  to  computed  affinities  fiom  the  Chem- 
score, Drug-score,  Flex-score,  G-score  and  PMF-score 
methods  respectively.  Although  an  improvement  was 
found  for  this  consensus  score  within  the  training  set  itself 
(R^  = 0.26),  its  predictive  potency  is  poor,  judging  by  no 
evident  correlation  within  the  16  molecule  testing  set. 

To  help  verify  the  physical  sensibility  of  our  model, 
we  have  mapped  out  the  residues  that  contribute 
significantly  to  the  scoring  function.  Coefficients  for  the 
20  most  important  residues  in  terms  of  electrostatic  and 
VDW  contributions  are  shown  in  Figure  1.  Negative 
coefficients  in  this  figure  can  be  interpreted  as  indicating 
ligand-residue  interactions  that  enhance  the  binding 
affinity,  while  positive  coefficients  correspond  to 
interactions  that  diminish  the  affinity.  In  terms  of 
absolute  magnitude,  Trp86,  Ile451,  Gly448,  Tyr449, 
Ser229  are  the  most  important  residues  in  the  active  site 
for  VDW  interactions.  Trp86  enhances  ligand  binding 
affinity  by  forming  n-n  interaction  with  ligand  aryl 
groups  (when  available),  while  the  other  residues  play  a 
largely  steric  role  by  defining  the  shape  of  the  gorge  base, 
thus  serving  to  discriminate  according  to  ligand  shape.  In 


22 


the  upper  gorge  area,  residues  Tyrl24,  Phe295,  Phe338, 
Phe297  are  responsible  for  providing  hydrophobic 
contacts.  The  ring  of  Tyr72  is  almost  perpendicular  to  the 
Trp286  ring  and  forms  a blocking  wall  to  prevent  the 
ligand  ring  moving  away  from  the  position  forming  k-k 
interaction  with  the  Trp286  ring.  Phe295,  Phe297, 
Val365  and  Glu292  form  another  wall  on  the  other  side  of 
the  PAS. 

Residues  Tyr449,  Glu450,  De451,  Alal27,  Serl28, 
Tyrl33,  Del  18  near  Try86  and  the  “oxyanion  hole” 
residues  Glyl21,  Glyl22  are  important  in  providing 
electrostatic  interactions  in  the  active  site.  Tyr337, 
Asp74,  Thr83  and  Asn87  are  the  primary  electrostatic 
contributors  in  the  gorge  area.  Gly342,  Lue76,  Glu285, 
Trp286,  His287  and  Gln291  likely  help  to  enhance  the 
activity  of  ligands  with  polar  groups  oriented  in  this  area, 
as  evidenced  by  reports  that  an  AChE  inhibitor  tethering 
in  the  position  of  His287  can  affect  the  binding  affinity  as 
much  as  14-fold.  Site-directed  mutagenesis  in  huAChE 
indicated  that  Asp74,  Tyr337,  Phe338,  Phe295,  Phe297, 
Tyrl33,  Glu450  can  affect  the  affinity  although  it  has 
been  difficult  to  determine  experimentally  whether  these 
residues  contribute  mainly  electrostatic  or  VDW 
interactions. 

The  docked  ligand  structures  generally  support  the 
above  analysis  regarding  the  identity  or  principle  residues. 
In  our  current  docking  calculations,  molecules  2-7  all 
share  similar  conformations  in  the  PAS.  The  modified 
groups  in  those  molecules  are  actually  exposed  in  the 
solvent  and  do  not  contribute  directly  to  the  ligand- 
receptor  interaction.  This  is  the  same  as  observed  in 
previous  studies**^.  However,  in  the  current  work, 
molecule  8 takes  a slightly  different  conformation  v«th  its 
morpholino  group  situating  in  the  half-buried  packet  by 
Trp286,  His287,  Ser293,  Glu292  and  Leu289.  The 
morpholino  Oxygen  has  a distance  of  3.65  A to  the 
backbone  N of  Glu292,  and  the  morpholino  Nitrogen  is 
3.78  A from  the  backbone  O of  Ser293.  Such  interactions 
could  slightly  pull  the  benzisoxazole  ring  away  from 
Trp286  ring,  leaving  only  a partial  K-n  interaction.  This 
particular  conformation  leads  to  an  affinity  increased  by 
more  than  10-folds  to  0.8nM  relative  to  molecules  2-7. 
This  particular  region  has  been  confirmed  in  an  x-ray 
structure  to  be  very  important  for  the  inhibitors  binding  to 
the  PAS  (periferal  anionic  site). 

As  a final  point  of  validation,  we  compared  the 
calculated  conformation  for  E2020  within  the  huAChE 
crystal  structure  relative  to  its  original  co-crystallized 
conformation  in  tcAChE.  Our  calculated  structure 
suggests  that  E2020  has  similar  but  not  identical  binding 
modes  in  tcAChE  and  huAChE.  In  the  active  site,  its 
benzyl  ring  forms  a n-n  interaction  with  the  indole  ring  of 
Trp86  in  huAChE  and  Trp84  in  tcAChE.  In  the  PAS,  the 
indanone  ring  of  E2020  forms  a tt-tt  interaction  with  the 
indole  ring  of  Trp286  in  huAChE  and  Trp279  in  tcAChE. 


The  charged  nitrogen  of  the  E2020  piperidine  ring 
undergoes  a cation-7i  interaction  with  the  phenyl  ring  of 
Phe330  in  tcAChE.  The  corresponding  residue  Tyr337  in 
huAChE  does  not  form  a similar  cation-7t  interaction  due 
to  the  steric  limitations  in  this  area  associated  with  an 
inauspicious  orientation  of  the  Tyr337  ring.  As  a result, 
the  nitrogen  of  the  piperidine  ring  of  the  E2020  instead 
interacts  with  the  hydroxyl  groups  of  the!Tyr337,  Tyrl24 
and  Serl25  within  distances  of  3.41,  3112  and  4.18  A 
relative  to  the  oxygen  atom  of  the  respective  hydroxyl 
groups.  To  probe  the  role  of  Tyr337,  we  examined  the 
potential  energy  curve  for  Tyr337  side  chain  rotation 
relative  to  the  other  residues  (energy  evaluation  according 
to  the  MMFF94s  force  field  with  appropriate  charges), 
but  found  there  to  be  only  one  minimum'  in  the  potential 
curve.  Closer  analysis  reveals  that  that  the  Tyr337  ring  is 
trapped  in  a local  pocket  formed  by  Phe338,  Tyr341, 
Trp439,  Trp449  and  His447.  In  previous  structural 
studies  and  molecular  dynamics  simulations,  it  has  been 
found  that  Phe330  in  tcAchE  can  adopt  a wide  range  of 
conformations  in  the  complex  structure  and  may  function 
as  a swinging  gate  that  structurally  couples  the  anion  sub- 
site of  the  active  site  and  the  PAS.  It  is  natural  to  expect 
similar  behavior  of  Tyr337  in  the  huAChE  as  compared 
with  the  analogous  Phe330  in  tcAChE.  Such  a gate  swing 
movement  of  the  Tyr337  ring  would  certainly  induce  a 
shift  in  attached  and  neighboring  residues,  thus  one 
function  of  this  PAS/active  site  coupling  could  be  to 
effect  proper  residue  alignment  within  the  anionic  subsite. 
The  fact  that  the  huAChE  crystal  structure  used  in  this 
study  to  train  our  scoring  function  did  |not  have  a co- 
crystallized active  site  inhibitor  (having  only  a PAS- 
bound  fasciculin  molecule)  and  thus  did  not  reflect  such  a 
conformational  shift  in  the  Tyr337  may  constitute  a subtle 
flaw  in  our  model.  Given  our  model's  fairly  strong 
predictive  capacity,  we  expect  the  flaw|  to  be  of  only 
minor  consequence,  however.  ] 

Competitive  Transport  Kinetics:!  The  Torpedo 
Califomica  AChE  receptor  cavity  entails  one  hydrophobic 
end  demarcated  by  Trp  84,  and  a polar  side  including  His 
440  and  Ser  200  that  initiates  the  hydrolysis  reaction. 
Molecules  such  as  ACh  and  many  of  the  OP  nerve  agents 
have  both  hydrophobic  and  polar  groups,  thus  one  can 
envision  a possible  binding  mechanism  such  that  the 
hydrophobic  end  of  the  ligand  is  attracted  to  the  gorge 
region  by  the  lipophilic  periferal  site,  migrates  to  the 
deeper  Trp  84  residue,  and  then  swings  its  polar 
electrophilic  group  down  toward  Ser  200.  The  tendency 
toward  such  behavior  can  be  measured  in  MD  simulations 
by  tracking  center  of  mass  movements  j throughout  the 
trajectory.  The  progression  of  a soman  molecule  has  been 
monitored  throughout  the  1 ns  simulation  (Ser  200  treated 
anionically).  This  initial  exterior  starting  position  leads  to 
rapid  (—30  ps)  admission  into  the  gorge,  at  which  point  its 
hydrophobic  1,2-dimethyl-propyl  group!  remains  fairly, 
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well  coupled  to  Trp  84.  A similar  simulation  of  DFP 
yielded  almost  identical  behavior,  strongly  suggesting  that 
relatively  small  molecules  such  as  ACh  and  most  OP 
nerve  agents  are  drawn  rapidly  into  the  cavity. 

Once  in  the  cavity,  however,  little  evidence  is  found 
for  strong  coupling  between  the  ligand’s  polar 
phosphonate  group  and  Ser  200.  This  is  indeed  case  for 
all  of  the  simulations  undertaken  with  the  deprotonated 
Ser  200.  This  is  further  reflected  in  studying  the  distances 
(Rco)  beween  the  ester  carboxyl  carbon  of  ACh  and  the 
side  chain  oxygen  of  Ser  200.  The  Rco  arising  from  the 
trajectory  are  summarized  according  to  a histogram  of 
ln[Frequency(Rco)]  versus  intervals  of  Rco  values.  This 
histogram  can  be  fitted  very  well  to  the  following 
relationship:  In(Rco)  = -4.077  Rco^  + 46.866  Rco 

-130.032.  This  can  be  used  to  predict  the  approximate 
frequency  whereby  Rco  closes  to  within  a plausible 
covalent  bonding  distance  (e.g.,  1.5  A):  approximately 
1.2x10“^”  instances/ns,  or  a rate  of  about  once  every 
2.7xl0‘^  years.  This  highly  unrealistic  finding  suggests 
that  Ser  200  deprotonation  must  not  be  a precursor  to  the 
Ser(0)  nucleophilic  attack  on  the  ester  (or  OP)  center:  an 
anionic  Ser(0)  likely  experiences  too  much  electrostatic 
repulsion  from  the  ester  (or  OP)  oxygens  to  permit 
nucleophilic  attach  on  the  ester  (OP)  center.  A neutral 
Ser  200  gives  a much  more  plausible  result.  In  this  case 
the  logarithmic  frequency  obeys  two  very  distinct  and 
separate  curves  of  which  the  one  on  the  right  can  be 
interpreted  as  a lipophilically  bound  conformer  attached 
primarily  to  Trp  84,  and  the  other  a structure 
characterized  by  both  lipophilic  Trp  84  and  polar  Ser  200 
+ His  440  binding.  The  latter  curve  permits  much  closer 
interaction  between  the  ligand’s  carbonyl  carbon 
electrophile  and  the  Ser(0)  nucleophile,  with  a 
logarithmic  frequency  described  by  In(Rco)  = -3.445 
Rco^  -0.334  Rco  + 3.955.  By  this  relationship,  the  same 
1.5  A separation  is  predicted  to  occur  at  a rate  of  0.1  / ns. 
Comparing  this  rate  with  AChE’s  known  overall  ACh 
processing  rate  of  ~10’  molecules  per  second  (i.e., 
0.001/ns),  and  recalling  the  rapidity  whereby  periferally 
bound  molecules  can  enter  into  the  gorge  (e.g..  Fig.  10), 
it  is  our  assessment  that  intragorge  diSlision  is  not  a rate 
determining  step  in  the  ACh/AChE  hydrolysis  reaction. 
Our  simulations  of  OP  compounds  indicate  that  they 
exhibit  diffusion  behavior  that  is  similar  to  that  of  ACh. 
We  find  that  sarin,  for  example,  exhibits  the  same  sort  of 
bimodal  logarithmic  frequency  plot  as  ACh. 
Interestingly,  although  it  spends  a significantly  greater 
amount  of  time  in  the  a bindable  conformation  at  a rate  of 
approximately  once  every  155  ps  - a value  that  is  much 
greater  than  that  of  ACh,  but  is  still  less  than  the  total 
time  required  for  most  OP  agents  to  phosphorylate  the 
AChE  enzyme  (in  the  ms  - minute  time  frame  for  ligand 
concentrations  equivalent  to  fatal  doses). 


Receptor  vibrations  that  enhance  ligand  transport 
behavior  have  been  observed  in  prior  simulations  of  the 
AChE  receptor.  Assuming  that  the  principle  function  of 
such  sympathetic  vibrations  is  to  eliminate  log-jam 
scenarios  that  would  otherwise  prevent  the  ligand  from  its 
appropriate  intra-receptor  transport,  we  have  sought  to 
characterize  such  phenomena  through  analysis  of  rapid 
jumps  in  the  root  mean  square  deviation  (RMSD)  of 
ligand  atom  positions  relative  to  the  starting  conformer. 
Ligand  RMSD  profiles  have  thus  been  plotted  and 
analyzed.  The  most  dramatic  shift  in  conformation  occurs 
for  VX  between  the  100  ps  and  200  ps  points  of  the 
simulation.  One  finds  that  the  VX  shift  (a  net  RMSD  of 
1.77  A)  corresponds  largely  to  a lurch  of  more  than  1.50 
A away  from  the  periferal  site  residue  Phe  33 1 and  closer 
to  the  active  site  Ser  200.  Despite  the  magnitude  of  this 
motion,  the  RMSD  between  the  corresponding 
conformations  of  Phe  330  + Phe  331  is  actually  very 
comparable  (1.74  A),  moving  substantially  out  of  the  way 
of  VX  in  its  approach  to  the  receptor.  While  this  shift  in 
these  two  residues  is  clearly  the  key  feature  of  any 
concerted  motion  in  the  protein,  the  rest  of  the  system 
does  also  shift  appreciably  (net  RMSD  = 1.20  A),  thus 
supporting  the  observation  of  Shen  et  al.^’'*'  who 
suggested  broad  scale  rearrangement  in  AChE  during 
ligand  admission  to  (or  expulsion  from)  the  gorge. 

Prediction  of  Adsorption  Kinetics, 
Thermodynamics,  and  Physical  Properties:  Much 

effort  has  been  devoted  to  painstaking  analysis  of  the  role 
played  by  the  components  of  the  catalytic  triad  and  the 
surrounding  residues  (individually  and  in  concert),  and 
their  effect  on  the  catalytic  power  of  the 

acetylcholinesterase  system.  Recent  proton  NMR 
studies’’^*  have  suggested  that  formation  of  a short  strong 
hydrogen  bond  (SSHB)  within  the  catalytic  triad  is  a key 
factor  in  the  efficiency  of  both  AChE  and  the  related 
enzyme  butyrylcholinesterase  (BChE).  These  studies 
have  demonstrated  a pronounced  downfield  shift 
(>15ppm)  in  the  IH  NMR  signal,  as  well  as  hydrogen 
bond  distances  <2.65  Angstroms.  We  have  expanded 
upon  oiu  previous  efforts  to  determine  the  exact  factors 
contributing  to  this  downfield  shift  upon  binding,  with 
surprising  results.  Table  1 contains  proton  NMR  shift  and 
residue  distances  obtained  from  our  quantum  calculations 
for  binding  (proton  and  aged  DEFP)  in  three  separate 
models  of  varying  size  and  complexity.  The  smallest 
model  contains  only  the  active  site  residue  sidechains 
(ACTside),  used  in  much  of  the  earlier  work  (ours  as  well 
as  others’)  on  this  system.  The  medium  sized  model 
contains  the  full  active  site  residues  as  well  as  the  full 
oxyanion  hole  residues  (ACT  + OX).  The  final,  largest 
model  expands  upon  the  medium  model  with  the  addition 
of  the  full  neighboring  residue  Glul99,  which  has  been 
postulated  to  play  a variety  of  roles  in  this  system.  A 
close  analysis  of  the  distances  within  the  active  site  (the 
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Glu327  - His  440  hydrogen  bond  distance)  brings  to  light 
very  interesting  details  of  this  reaction.  The  experimental 
data  suggests  a FURTHER  downfield  shift  in  the  proton 
NMR  upon  binding,  specifically  from  14.4  ppm  for  bare 
acetylcholinesterase  to  15.5  and  16.6  for  paraoxon  (an 
analogue  of  DEFP).  This  coincides  with  an  actual 
contraction  of  the  hydrogen  bonds  across  the  active  site 
residues.  The  exact  numerical  values  obtained  from  these 
models  do  not  invite  exact  comparison  with  experimental 
data,  as  NMR  data  is  notoriously  difficult  to  obtain  to 
quantitative  accuracy.  However,  the  values  obtained  are 
all  within  reasonable  accuracy  for  the  method  and  basis 
set  size  used,  and  shows  extremely  meaningful  trends.  To 
whit,  the  only  model  which  displays  tightening  of  the 
active  site  and  a correct  downfield  shift  of  the  GluH 
signal  relative  to  the  bare  system  is  the  largest  model 
containing  the  Glul99,  which  in  this  run  has  formed  a 
CHO  hydrogen  bond  to  the  histidine  ring  (see  Figure  2). 
This  conclusively  demonstrates  the  importance  of  this 
residue  in  the  binding  process,  as  postulated 
experimentally. 


Table  1. 


System 

Substrate 

Glu327- 

His440 

distance 

GluH 

shift 

HisH 

shift 

ACTside 

X 

2.614 

18.972 

8.039 

ACTside 

H 

2.712 

14.453 

11.207 

ACT+OX 

X 

2.584 

19.717 

12.722 

ACT+OX 

Aged  DEFP 

2.651 

16.938 

14.391 

ACT+OX 

H 

2.693 

14.952 

11.095 

ACT+OX+GLU1 

99 

X 

2.623 

18.854 

13.022 

ACT+OX+GLU1 

99 

Aged  DEFP 

2.549 

21.857 

15.601 
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Residues 


a)  Residues  importance  in  VDW  interactions 


Figure  2.  ACT-«-OX+Glu199  with  bound  agent 
stimulant,  demonstrating  CHO  hydrogen  bonding  of 
Glu199  to  HIS440 


Residues 


b)  Residues  importance  in  electrostatic  interactions 

Figure  1.  a)  Coefficients  of  the  20  most  important 
residues  for  electrostatic  interactions,  b)  Coefficients 
of  the  20  most  important  residues  for  VDW 
interactions. 
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